################################################################################
##############17.Robustness Check with WIID#####################################
################################################################################

#####Load data from Data_Generation_CPS#####

pols <- readRDS('pols_bjps')


#####Load WIID####

library(readxl)
wiid <- read_excel("WIID_17Dec2019 (1).xlsx") %>%
  mutate(country = ifelse(country == 'United Kingdom', 'Great Britain',
                               country),
         old = 0) %>%
  rename(country.name = country,
         electionyear = year) %>%
  filter(source == 'Eurostat' | source == 'Luxembourg Income Study' | source == 'OECD' | source == 'World Bank') %>%
  filter( resource == 'Income (net)' & scale == 'Equivalized' & scale_detailed == 'OECD-modified' & quality == 'High' & areacovr_detailed == 'All')%>%
  dplyr::select(country.name, electionyear, gini_reported, resource, resource_detailed, scale, scale_detailed, quality)

pols <- pols %>%
  mutate(old = 1)

pols_wiid <- left_join(pols, wiid, by = c('country.name', 'electionyear')) %>%
  mutate(scale_gini = scale(gini_reported)) %>%
  filter(old ==1)

library(plm)
pols_wiid <- pdata.frame(pols_wiid, index = c('country.name', 'electionyear'))

cor.test(pols_wiid$gini_disp, pols_wiid$gini_reported)

wiid1 <- plm(scale_econsdw ~ scale_gini + 
               scale_ggini_inc1 + 
               scale_socioeconomicsal +
               scale_gini : scale_ggini_inc1 + 
               scale_gini:scale_socioeconomicsal +
               scale_ggini_inc1: scale_socioeconomicsal +
               scale_ggini_inc1:scale_socioeconomicsal : scale_gini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols_wiid, 
             na.action = na.omit, 
             model = "within")

coeftest(wiid1, vcov = function(x) vcovHC(x))

wiid2 <- plm(scale_econsdw ~ scale_gini + 
               scale_ggini_inc2 + 
               scale_socioeconomicsal +
               scale_gini : scale_ggini_inc2 + 
               scale_gini*scale_socioeconomicsal +
               scale_ggini_inc2: scale_socioeconomicsal +
               scale_ggini_inc2:scale_socioeconomicsal : scale_gini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols_wiid, 
             na.action = na.omit, 
             model = "within")

coeftest(wiid2, vcov = vcovHC)

wiid3 <- plm(scale_econsdw ~ scale_gini + 
               scale_ggini_inc3 + 
               scale_socioeconomicsal +
               scale_gini : scale_ggini_inc3 + 
               scale_gini*scale_socioeconomicsal +
               scale_ggini_inc3: scale_socioeconomicsal +
               scale_ggini_inc3:scale_socioeconomicsal : scale_gini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols_wiid, 
             na.action = na.omit, 
             model = "within")

coeftest(wiid3, vcov = vcovHC)

wiid4 <- plm(scale_econsdw ~ scale_gini + 
               scale_ggini_inc4 + 
               scale_socioeconomicsal +
               scale_gini : scale_ggini_inc4 + 
               scale_gini*scale_socioeconomicsal +
               scale_ggini_inc4: scale_socioeconomicsal +
               scale_ggini_inc4:scale_socioeconomicsal : scale_gini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols_wiid, 
             na.action = na.omit, 
             model = "within")

coeftest(wiid4, vcov = vcovHC)

wiid5 <- plm(scale_econsdw ~ scale_gini + 
               scale_ggini_inc5 + 
               scale_socioeconomicsal +
               scale_gini : scale_ggini_inc5 + 
               scale_gini*scale_socioeconomicsal +
               scale_ggini_inc5: scale_socioeconomicsal +
               scale_ggini_inc5:scale_socioeconomicsal : scale_gini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols_wiid, 
             na.action = na.omit, 
             model = "within")

coeftest(wiid5, vcov = vcovHC)

##### Table 1: Model 4 Output#####

fit1c <- c(summary(wiid1)$r.squared, summary(wiid1)$fstatistic$statistic, summary(wiid1)$fstatistic$p.value)
fit2c <- c(summary(wiid2)$r.squared, summary(wiid2)$fstatistic$statistic, summary(wiid2)$fstatistic$p.value)
fit3c <- c(summary(wiid3)$r.squared, summary(wiid3)$fstatistic$statistic, summary(wiid3)$fstatistic$p.value)
fit4c <- c(summary(wiid4)$r.squared, summary(wiid4)$fstatistic$statistic, summary(wiid4)$fstatistic$p.value)
fit5c <- c(summary(wiid5)$r.squared, summary(wiid5)$fstatistic$statistic, summary(wiid5)$fstatistic$p.value)


fit.statsc <- rbind(fit1c, fit2c, fit3c, fit4c, fit5c)

mean(fit.statsc[,1]) #R-Squared
mean(fit.statsc[,2]) #Adjusted R-Squred
mean(fit.statsc[,3]) #F-Statistic
mean(fit.statsc[,4]) #p-value of F-Statistic

fu1 <- as.data.frame(summary(wiid1, vcov = vcovHC)$coefficients) 

fu2 <- as.data.frame(summary(wiid2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(wiid3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(wiid4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(wiid5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

#dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
#  geom_point() + 
#  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
#  geom_linerange(aes(ymin = lwr, ymax = up))+
#  coord_flip() +
#  labs(title = "Point Estimates of Model Regressing Predictors #on Economic Polarization (n = 83)", y = "Point Estimates (95% #Confidence Interval)", x = "Predictor")+
#  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

pols_rug <- pols_wiid %>%
  dplyr::select(scale_econsdw,scale_gini, 
                scale_ggini_inc1,
                scale_socioeconomicsal,
                scale_galtansdw ,
                scale_enep ,
                scale_unemp) %>%
  na.omit

test <- as.data.frame(mvrnorm(1000, mu =coef(wiid1), Sigma =vcovHC(wiid1)))

#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols_wiid$scale_ggini_inc1, na.rm = T), to = max(pols_wiid$scale_ggini_inc1, na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*-1 + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sort <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

lwr.plot.v <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = -1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

#These steps are repeated for each of the panels in each figure.#

#####Middle Panel####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*0 + test[j,10]*marg[i,1]*0
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

middle.plot.v <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 0", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

####Right Plot#####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*marg[i,1] + test[j,8]*1 + test[j,10]*marg[i,1]*1
  }
}

sort <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sort = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

righ.plot.v <- ggplot(marg.plot.dat, aes(x = sort, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Economic Salience = 1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_ggini_inc1), inherit.aes = F)

wiid_sortx <- grid.arrange(lwr.plot.v, middle.plot.v, righ.plot.v, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Sorting by Income")

#ggsave("figure_A5.pdf", plot = wiid_sortx, width = 11, height = 5, units = 'in')

#####Figure 3: Salience on x-axis, sorting held constant#####


#vector of sorting from observed min to observed max#
nums <- seq(from = min(pols_wiid$scale_socioeconomicsal, na.rm = T), to = max(pols_wiid$scale_socioeconomicsal, na.rm = T), length.out = 1000)

#create data frame for storing output#
marg <- data.frame(x = nums)

#Simulate marginal effect based on draws from the model#
for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*-1 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*-1
  }
}

#Vectors to be filled#
sal <- NA
pe <- NA
lb <- NA
ub <- NA


#Make CI based on empirical quantiles

for(i in 1:1000){
  sal[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

#Combine into a singe dataset

marg.plot.dat <- data.frame(sal = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

#create plot. 

lwr.plot.v <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = -1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)

#These steps are repeated for each of the panels in each figure.#

#####Middle Panel####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*0 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*0
  }
}

sal <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sort[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(sal = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

middle.plot.v <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = 0", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)


####Right Plot#####

for(i in 1:1000){
  for(j in 1:1000){
    marg[i,j+1] <- test[j,1] + test[j,7]*1 + test[j,8]*marg[i,1] + test[j,10]*marg[i,1]*1
  }
}

sal <- NA
pe <- NA
lb <- NA
ub <- NA



for(i in 1:1000){
  sal[i] <- marg$x[i]
  pe[i] <- unlist(quantile(marg[i, 2:1001], probs = .5))
  lb[i] <- unlist(quantile(marg[i, 2:1001], probs = .025))
  ub[i] <- unlist(quantile(marg[i, 2:1001], probs = .975))
}

marg.plot.dat <- data.frame(al = nums,
                            pe = unlist(pe),
                            lb = unlist(lb),
                            ub = unlist(ub))

righ.plot.v <- ggplot(marg.plot.dat, aes(x = sal, y = pe)) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = 'grey69', alpha = .5) +
  geom_line() +
  ylim(c(-1.5, 3.5))+
  geom_hline(aes(yintercept = 0), linetype = 'longdash', color = 'red') +
  labs(title = "Sorting = 1", x = "", y = "") +
  theme_bw() +
  geom_rug(data = pols_rug, aes(x = scale_socioeconomicsal), inherit.aes = F)


wiid_salx <- grid.arrange(lwr.plot.v, middle.plot.v, righ.plot.v, nrow = 1, top = "", left = "Marginal Effect of Disposable Income Inequality", bottom = "Salience of Economic Issues")

#ggsave("figure_A6.pdf", plot = wiid_salx, width = 11, height = 5, units = 'in')

